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ABSTRACT 

We study the escape of cosmic-ray (CR) protons accelerated at a supernova remnant 
(SNR) by numerically solving a diffusion-convection equation from the vicinity of the 
shock front to the region far away from the front. We consider the amplifications of 
Alfven waves generated by the escaping CR particles and their effects on CR escape 
into interstellar medium (ISM). We find that the amplification of the waves signif- 
icantly delays the escape of the particles even far away from the shock front (on a 
scale of the SNR). This means that the energy spectrum of CR particles measured 
through 7-ray observations at molecular clouds around SNRs is seriously affected by 
the particle scattering by the waves. 
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1 INTRODUCTION 

Cosmic rays (CRs) are relativistic charged particles and 
most of them are believed to be accelerated at the shock 
front of supernova remnants (SNRs). In fact, TeV 7-ray 
emissions have b een detected with H.E. S .S. fr om the shell 
of young SNRs (|Aharonian et al ] l20o4 I2005I ). Moreover, 
GeV 7-ray emissions have been detected with Fermi and 
AGILE from middle-aged SNRs surrounded by molecular 
clouds. They are likely to be generated via a hadronic pro- 
cess, i.e. inelastic collisions bet ween CR protons accelera ted 
at SNRs and ambi e nt protons (lAbdo et al.ll2009l . I2010allbll3 : 
iTavani et all I2OI0I : ICiuhani et al.ll2010l ). TeV 7-ray emis- 
sions have also been d iscovered from some of these clouds 
lAharonian et alll2008l ). The energy spectrum of CR pro- 
tons is not described by a single power-law. Thus, it is most 
likely to be affected by an e scape process from SNR shocks 
(but see also lUchivama et~al.,,2010, V 

The 7-ray observations of molecular clouds around 
SNRs have given us information not only on the 
CR particle acceleration in the SNRs but also on 
the escape of the particles from the SNRs (e.g. 
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their escape into the surrounding interstellar medium 
(ISM) is important to compare the 7-ray observations of 
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molecular clouds around the SNRs with theoretical models 
of pa rticle acceleration. Recently, Lee. Kamae. fc EllisonI 
Hooi) and ICabici. Aharonian. fc Casanoval (|2009l 'l theo- 
retically studied the escape into the ISM and obtained an 
energy spectrum of the escaped particles. In their studies, 
they assumed that the diffusion coefficient for the particle 
diffusion outside the SN R is the one in th e general region 
in the Galaxy. However, iFuiita et al.l (|2009l ) suggested that 
the observations seem to show that the diffusion time-scale 
of CRs around SNRs is much longer than that in the 
general region in the Galaxy. They estimated that the 
diffusion coefficient must be less than ~ 1 % of that in the 
general region to be consistent with the observations of 
the SNR W 28 and another possible SNR in Westerlund 2 
(see also [Torres. Rodriguez Marrero. fc de Cea Del Pozc 
I2OO3). Motived by these results, iFuiita. Ohira. fc Takahara 
(2010|) theoretically studied the diffusion of particles in 
ISM after they left the shock neighbourhood. They used 
Monte-Carlo simulations to follow a particle motion, while 
at the same time they calculated the evolution of the SNR. 
They considered the excitation of Alfven waves through the 
streaming instability caused by CRs, and the interactions 
between the waves and the particles. Since the growth of 
the waves changes the diffusion coefficient, they took the 
evolution of the diffusion coefficient into consideration. 
They found that the particles actually excite Alfven waves 
around the SNR on a spatial scale of the SNR itself if the 
ISM is highly ionised. Thus, even if the particles can leave 
the shock neighbourhood, scattering by the waves prevents 
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them from moving further away from the SNR. This means 
that the particles cannot virtuahy escape from the SNR 
until a fairly late stage of the SNR evolution. 



However, 



the 



model 



developed 



by 



iFuiita. Ohira. fc Takaharal (|2010l) may still be rather 
simplified. In particular, they separately treated the accel- 
eration of particles at the shock front and their propagation 
into the surrounding ISM. Since both of the processes 
follow the same transport equation, they should be treated 
seamlessly. In this paper, we present the results of our 
numerical simulations that deal with both the particle 
acceleration and the diffusion into the ISM at the same 
time. In these simulations, we take account of the streaming 
instability developed by the flux of CRs in order to study 
how the growth of the Alfven waves affects the diffusion of 
the CRs. We emphasise that we focus on the propagation 
of CRs far away from the shock front {r — Rs ~ Rs, 
where r is the distance from the SNR centre and -Rs is 
the shock radius), in contrast with previous studies that 
treat the escape around the shock f ront (r — i?s <^ Rs; e.g. 
IVladimirov, EUison. fc Bvkovl |2006| ). As far as we know, 
this is the first time to solve a transport equation from the 
shock vicinity to the region far away from the shock front. 



2 MODELS 
2.1 Equations 

The diffusion of CR particles obeys the transport equation 
(diffusion-convection equation): 



(1) 



where f{r,p,t) is the distribution function of the parti- 
cles, p is the momentum, k is the diffusion coefficient, w 
is the hydrodynamic velocity of the background gas, and 
Q = QoS{r — Rs) is the source term for the particles, which 
are injected at the shock front (r = Rs) at a momentum of 
Pinj. The coefficient Qo is given by 
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where pi and u\ are the gas density and the gas velocity 
relative to the shock just upstream the shock front, respec- 
tively. In § O we set pinj ~ 2 mcs2, where Cs2 is the sound 
velocity behind the shock front. The fraction of gas parti- 
cles that go into the acceleration process at the shock front 
is e = 10"". 

The equations for the background gas are 
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where p, 7g, and Pg are the density, specific heat capacity 
ratio, and pressure of the gas. The CR pressure is given by 
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where Pmin is the minimum momentum of injected CR par- 
ticles, m is the mass of the particles, and c is the velocity of 
light. We consider only protons and neglect electrons. 



2.2 Diffusion Coefficient and Initial Conditions 

Since we treat the propagation of CRs away from the 
shock front, we cannot assume a plane geometry. There- 
fore, we assume spherical symmetry and solve equations 
([T]) -([51) mainly based on the numeric al methods developed 
bv lBerezhko. Yelshin. fc Ksenofonto'^ ([l994). We do not fix 
the diffusion coefficient k at the value corresponding to 
the Bohm diffusion, which is diffe r ent fr om the studies by 
iBerezhko. Yelshin. fc Ksenofontovl (| 19941 '). Instead, we fol- 
low the evolution of k in this study, considering the ampli- 
fication of Alfven waves by CR particles. Since the waves 
scatter the particles, they affect the diffusion coefficient in 
equation IT}. The growth of the waves is given by 
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where xp{t,r,p) is the energy density of Alfven waves per 
unit logarithmic bandwidth (which are resonant with parti- 
cles of momentum p) relative to the ambient magnetic en- 
ergy dens ity Uyi, and va is the Alfven velocity (|Skillinell 19751 : 
[Bell 19781 ). We do not consider the damping of the waves for 
simplicity. The diffusion coefficient is 
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where v is the velocity of the particle, e is the ele- 
mentary charge. Bo is the unperturbed magnetic field, 
and Po is the density of unperturbed ISM. Following 
Ber ezhko. Yelshin. fc Ksenofo ntov (1994), we included the 
infiuence of gas compression on the magnetic field {B = 
Bop/ Po) and particle scattering by assuming that k depends 
on the gas density. 

The diffusion coefficient corresponding to the Bohm 
diffusion is given by kbo = Pbc/3, where pb = pc/{eBo) 
is the gyroradius of particles. We calculate the evolu- 
tion of the gas and CR particle s usin g the methods of 
IBerezhko. Yelshin. fc Ksenofontovl (|l994l ) for r < ri,{p) = 
4kbo/V's -I- Rs and p/{mc) > 10~^ *, where Vs is the shock 
velocity. We call this region the "shock region". Since we 
are also interested in the CR diffusion outside this region, 
we added equal interval meshes for rb{p) ^ r < rout(p) = 
500 hbo/Vs + rb. We call this surrounding region the "es- 
cape region". For this region, we calculate the diffusion of 
the CR particles only with p/(mc) > 10^ *. For this energy 
range, the escape region is sufficiently apart from the shock 
front, and the feedback of the CR particles on hydrodynam- 
ics is small there. Thus, we ignore the feedback in this region 
(Viy ~ 0) and solve equation ((T} without considering the 
change of / in the direction p (the third term in the right 
hand of equation [1]), which makes calculations much easier. 
The term wVf in equation ((TJ can also be ignored in this 
region, although we do not omit the term in the calculations. 
On the other hand, we do not solve the diffusion of parti- 
cles with p/(mc) < 10^ '* in the escape region, because Vto 
cannot be ignored even in the region. Since the feedback of 
the low-energy particles on hydrodynamics is limited in the 
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vicinity of the shock, they do not afTect the following results 
at r — Rs '-^ Rs ■ We set / = at r = rout • 

The initial condition of ip is given as follows. The diffu- 
sion coefficient of the background fSM is given by 



«ISM 



10^** cm^ s^ 



VfO GeV/ 



Bo 
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where E is the particle energy 

ijGabici. Aharonian. fc Casanoval l2009t ). We relate kism to 
the energy density of Alfven waves ('(/'ism): 



KiSM = 



pvc 



(10) 



If we set ip — '(/'ISM ai t — tinj, where tinj is the time 
when we start injecting CRs, ip does not grow enough and 
particles are not accelerated. As has been well known, for 
diffusive shock acceleration to work, ^ must be as l arge 
as that for Bohm diffus ion near the shock front (e.g. iBelll 
1 19781 : iLucek fc Belli [20001 . Thus, in this study, we simply 
assume that the diffusion at the shock front can be rep- 
resented by Bohm diffusion {k = kboPo/p)- One can find 
that the energy density of Alfven waves corresponding to 
kbo is tl^B ~ 4v/{t7c) ~ 4/7r. On the other hand, for 
r — i?s 3> hbo/Vs, the coefficient should be the one at the 
interstellar space (kq = kism). As the initial condition of tp 
at t = tinj , we interpolate the two tps: 



i^ii^jP) = (V'b - ^ism) exp 

+^ISM , 



(r - R,{Uni))Vs{tinj) 
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where is a parameter. Although the value of cannot be 
specified, it would not be much larger than one. Therefore 
we assume = 5 in the following simulations. The larger 
tti gives the smaller diffusion coefficient at a given r and 
results in less escape of the particles. For values of as 
small as one, it is increasingly difficult to obtain stable and 
precise results, because / changes rapidly at the boundary 
between the shock and escape region, where resolution is 
not high enough. We consider three models: (A) the growth 
of is considered, (B) the diffusion coefficient is fixed at 
K ~ kboPo/p regardless of time and position, and (C) the 
wave energy density does not change and is fixed at tp — ^i, 
replacing tinj with t in equation pip . The models B and C 
are calculated for comparison. In Model A, we set a lower 
limit of t/) as t/i ^ i/'i, replacing tinj with t in equation pip . 
Moreover, we set a upper limit of '0 as ■(/) 5C i/)b. 

Since we do not consider the neutral damping of the 
waves, we choose the ISM with a relatively high tempera- 
ture. The density and sound velocity of the background ISM 
is po = 7.0 X 10~^^gcm~^ and Cs = 154kms~^, respectively. 
A supernova explodes at t = and r — with an energy of 
10^^ erg. The background magnetic field is Bo ~ ^SttUm ~ 
3 /xG. We inject CR particles and calculate the growth of 
Alfven waves (equation [7]) for t > tinj — 0.55 to, where to 
is the time when the free expansion phase of the SNR ends 
and the Sedov phase begins. For the parameters we adopted, 
to = 4.7 X 10^ yr. 



3 RESULTS AND DISCUSSION 

In Figure [T^, we show the profiles of V at t = 10 to for 
Models A and C. At this time, the radius of the SNR 
is Rs = 73 pc, and the compression ratio of the gas is 
P^/po ~ 3.6, where p2 is the gas density behind the shock, 
for Models A, B, and C. The left ends of the curves corre- 
spond to ipB, and the right ends correspond to ^ism- The 
difference between Model A and C indicates that CR par- 
ticles can significantly amplify Alfven waves. In particular, 
the amplification occurs even in the region far away from the 
shock front (say, at r ~ 27?^) for particles with pc = 20 TeV. 
The profiles are given by the balance between the growth 
through streaming instability and the fiow of \ow-ip gas from 
upstream on the shock coordinate. The amplification of the 
waves prohibits the diffusion of particles. In Figure [ijj, we 
show the profiles of / at 10 to- While the particles are con- 
fined around the shock front in Model B, they can extend 
to the region where ip ~ '(/'ism in Model C (Figure [1} . The 
results when the growth of ip is considered (Model A) are 
located between those two models, although the profiles of 
ip and / are somewhat different between pc= 1 and 20 TeV 
(Figure [T]) . The overall features of tp and / do not change 
for t > to. 

The results for Model A can be explained as follows. 
Figure [2] shows the evolution of ip and / for Model A for 
tinj < t ^ to, normalised by Rs{t) in the radial direction. 
The radius of the SNR at t = to is Rs = 23 pc. For pc — 
1 TeV, tp increases in the "instep" or the region where the 
gradient of ip is small (0.02 < (r - Rs)/Rs < 2). This is 
because of the leak of particles into the \ow-tp region {tp ~ 
'(/'ism). On the other hand, for pc = 20 TeV, tp increases in 
the "shin" or the region where the gradient of tp is large 
(0-1 < {r - Rs)/Rs < 0.6). The leak of particles into the 
\ow-tp region {tp ~ '(/'ism) is less significant for pc = 20 TeV 
than for pc = 1 TeV (Figure [2]d) . This is because i/'ism at 
pc=l TeV is smaller than that at pc = 20 TeV (Figure [SJi 
and equations 1^ and llOjl . For smaller ?/'ism, particles tend 
to go father away from the shock front, if the distance is 
represented in the units of kbo/Vs- Thus, for pc — 1 TeV, 
particles escaping into the instep make the gradient of / in 
the instep, while for pc — 20 TeV most particles are confined 
in the shin and make the gradient of / in the shin. The 
gradient of / increases tp (equation [7]). 

CR particles that escape from an SNR illuminate molec- 
ular clouds around the SNR and could be observed in 7-rays 
through pion production. Figure |3] shows the spectra of par- 
ticles at r = 1.1 i?s and 2 .Rj, at t = 10 to. In Model B, only 
particles with high energies can reach r ~ 1.1 _Rs, because 
the diffusion coefficient for low-energy particles is small. On 
the other hand, in Models A and C, even low-energy parti- 
cles can reach this radius, and the trend is more prominent 
at larger radii (r ~ 2 Rs). Comparison between Models A 
and C shows that the effect of the wave growth cannot be 
ignored in the region away from the shock front {r = 2 Rs); 
the energy density of particles at p/mc 10*'^ in Model A 
is smaller than that in Model C by more than a factor of 
10. In Figure [3] the broad hump of p**/ at p/mc ~ lO'''^ at 
r ^ 2 Rs for Models A and C reflects the escape of parti- 
cles into the instep region oi tp ^ tpism shown in Figure [T] 



for pc = 1 TeV. The dip at p/mc 



lO'' 



at r ~ 2 Rs for 



Model A is made because the particles hardly reach that 
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{r-Rs)IRs 




{r - RsVRs 

Figure 1. Profiles of wave energy density at t = 10 to • Thin 
lines correspond to waves interacting with particles with pc = 
1 TeV and thick lines correspond to those with pc = 20 TeV. 
Solid lines are for Model A and dashed lines are for Model C. (b) 
Same as (a) but for CR distribution /. 

radius because of the growth of tp in the shin (pc = 20 TeV 
in Figure [T]). The peaks at p/mc ~ lO^'^ for Models A, B 
and C mean that the region of -i/j ~ tpB and the particles 
confined in that region reach r ~ 2 i?s and those particles 
are freely escaping in our treatment effectively. The decrease 
of spectra for Models A and C at r = 2 for p/mc < 10^'^ 
suggests that the spectra are affected by the outer bound- 
ary condition of / = 0, because rc,ut{p) is relatively close to 
2Rs. If rout(p) is set at a larger radius, we expect that more 
particles reach r > 2 i?s • 

At t = 10 to, the total energies of CRs for r > 2 Ra 
are 4.0 x W*, 2.3 x 10''^ and 3.0 x lO*'* erg for Models A, 
B and C, respectively. These results show that the growth 
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{r - RsVRs 

Figure 2. (a) Profiles of wave energy density tp for Model A at 
t = 0.55 to, 0.60 to, 0.65 to, 0.70 to, and 1.0 to (from bottom to 
top), (b) Same as (a) but for CR distribution / at t = 0.60 to, 
0.65 to, 0.70 to, and 1.0 to (from bottom to top). 



of the waves can significantly change the CR spectrum and 
density around a middle-aged SNR, which should affect the 
7-ray emission originated from the CRs. In this study, we did 
not treat particles with pc ~ GeV. However, above discus- 
sion suggests that many of them can escape into the instep 
region, because -i/jism for these energies are fairly small. It 
is to be noted that the results shown here do not involve 
the effect of wave damping. For example, if the ISM is not 
fully ionised, the Alfven waves may da mp through collisions 
between charged and neutral particles (|Kulsrud fc Cesarskvl 
ll97ll : IO'C Drurv. Duffy, fc Kirklll99d ). If this happens, the 
spectrum of the escaped CR particles wou ld be significantly 
changed l|Fuiita. Qhira. fc Takaharall2010l '). 
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Figure 3. Spectra of CRs at r = 1.1 Rs and 2 at t = 10 to- 
Solid, dotted, and daslied lines are for Models A, B, and C, re- 
spectively. 



O'C Drury L., Duffy P., Kirk J. G., 1996, A&A, 309, 1002 
Fujita Y., Ohira Y., Takahara F., 2010, ApJ, 712, L153 
Fujita Y., Ohira Y., Tanaka S. J., Takahara F., 2009, ApJ, 
707, LI 79 

Gabici S., Aharonian F. A., Casanova S., 2009, MNRAS, 
396, 1629 

GiuHani A., et al., 2010, A&A, 516, LU 
Kawanaka N., loka K., Ohira Y., Kashiyama K., 2011, ApJ, 
729, 93 

Kulsrud R. M., Cesarsky C. J., 1971, ApL, 8, 189 
Lee S.-H., Kamae T., EUison D. C., 2008, ApJ, 686, 325 
Lucek S. G., Bell A. R., 2000, MNRAS, 314, 65 
Ohira Y., Murase K., Yamazaki R., 2010, A&A, 513, A17 
Ohira Y., Murase K., Yamazaki R., 2011, MNRAS, 410, 
1577 

Ptuskin V. S., Zirakashvili V. N., 2003, A&A, 403, 1 

SkiUing J., 1975, MNRAS, 173, 255 

Tavani M., et al., 2010, ApJ, 710, L151 

Torres D. F., Rodriguez Marrero A. Y., de Cea Del Pozo 

E., 2008, MNRAS, 387, L59 
Uchiyama Y., Blandford R. D., Funk S., Tajima H., Tanaka 

T., 2010, ApJ, 723, L122 
Vladimirov A., EUison D. C., Bykov A., 2006, ApJ, 652, 

1246 



4 CONCLUSION 

We have investigated the escape of CR protons accelerated 
at an SNR and their diffusion in the surrounding ISM. We 
solved a diffusion-convection equation from the vicinity of 
the shock front to the region far away from the front. We also 
considered the amplifications of Alfven waves generated by 
the escaping CR particles, which affects the diffusion of the 
particles into the ambient ISM. We found that the amplifica- 
tion of the waves reduces the diffusion coefficient on a scale 
of the SNR and the escape of the particles is significantly 
delayed. Our results suggests that the particle scattering by 
the waves should have infiuence on the energy spectrum of 
CR particles observed at molecular clouds around SNRs in 
the 7-ray band. 
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